void Calculate_LU(std::vector<std::vector<double> > &_L, 
                  std::vector<std::vector<double> > &_U, 
                  std::vector<double> &_b);
void Calculate_Ly(std::vector<std::vector<double> > &_L, 
                  std::vector<double> &_b);
void Calculate_Ux(std::vector<std::vector<double> > &_U, 
                  std::vector<double> &_y);



void Calculate_LU(std::vector<std::vector<double> > &_L, std::vector<std::vector<double> > &_U, std::vector<double> &_b)
{
    Calculate_Ly(_L, _b);
    Calculate_Ux(_U, _b);
}

void Calculate_Ly(std::vector<std::vector<double> > &_L, std::vector<double> &_b)
{
    int i, j;
    int n = _L.size();
    for(j = 0; j <= n - 2; ++j)
    {
        _b[j] = _b[j] / _L[j][j];
        for(i = j + 1; i <= n - 1; ++i)
        {
            _b[i] = _b[i] - _b[j] * _L[i][j];
        }
    }
    _b[n - 1] = _b[n - 1] / _L[n - 1][n - 1];
}

void Calculate_Ux(std::vector<std::vector<double> > &_U, std::vector<double> &_y)
{
    int i, j, k;
    int n = _U.size();
    for(j = n - 1; j >= 1; --j)
    {
        _y[j] = _y[j] / _U[j][j];
        for(i = 0; i <= j - 1; ++i)
        {
            _y[i] = _y[i] - _y[j] * _U[i][j];
        }
    }
    _y[0] = _y[0] / _U[0][0];
}   